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We present an extensive analysis of systematic deviations 
in Wolff cluster simulations of the critical Ising model, using 
random numbers generated by binary shift registers. We in- 
vestigate how these deviations depend on the lattice size, the 
shift-register length, and the number of bits correlated by the 
production rule. They appear to satisfy scaling relations. 



The main advantage of cluster Monte Carlo algorithms 
is that they suppress critical slowing down |ffl,||. For 
this reason, cluster algorithms are being explored ex- 
tensively ||. This has even led to the construction of 
special-purpose processors using the Wolff cluster algo- 
rithm @§. 

The problem of generating random numbers of suffi- 
cient quality is known to be complicated since the first 
computer experiments [|| . Many of the widely used algo- 
rithms are of the shift-register (SR) type (7j. These are 
extremely fast g|, can be implemented simply in hard- 
ware and produce 'good random numbers' with an 
extremely long perio d fm . 

Ferrenberg et al. found that the combination of 
the two most efficient algorithms (the Wolff cluster algo- 
rithm and the shift-register random-number generator) 
produced large systematic deviations for the 2D Ising 
model on a 16 x 16 lattice (see also [^2| ). Also random- 
walk algorithms appeared to be sensitive to effects due 
to the random- number generator |l3j| . 

Remarkably, we did not find visible deviations in simu- 
lations [^|Jl^1 performed on the special-purpose processor 
with the Wolff algorithm and a Kirkpatrick-Stoll random- 
number generator for lattices larger than 256 x 256. 

Motivated by this paradoxical situation, we made an 
extensive analysis of this problem using SGI workstations 
at the Delft University and a DEC AXP 4000/620 server 
at the Landau Institute. A total of about two thousands 
hours of CPU time was spent. 

We find several interesting facts. First, the maximum 
deviations occur at lattice sizes for which average Wolff 
cluster size coincides with the length p of the SR. 

Second, the deviations obey scaling laws with respect 
to p: they can be collapsed on a single curve. This opens 
the possibility to predict the magnitude of the systematic 
errors in a given quantity, depending on the lattice size, 
the shift-register length and, to some extent, also on the 
number of terms in production rule. 

Third, the deviations change sign when we invert the 



range of the random number: x — ► 1 — x. This provides 
a simple test, in two runs only, for the presence of sys- 
tematic errors. 

Finally, we introduce a simple ID random-walker 
model explaining how the correlations in the SR lead to 
a bias in Monte Carlo results. 

As a first step in understanding the results, it is natural 
to compare the length scales associated with the Monte 
Carlo process and the random generator. The first char- 
acteristic length is the mean Wolff cluster size (c). The 
second characteristic length is the size p of the shift reg- 
ister. The production rule 



(1) 



where ® is the 'eXclusive OR' operation, leads to three- 
bit correlations over a length p. So, it not surprising that 
the largest deviations occur at the lattice size L max for 
which these two lengths coincide. Since the mean Wolff 
cluster size behaves [|| as the magnetic susceptibility \, 
we expect at criticality that 
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where 7 and v are the susceptibility and correlation 
length exponents respectively. 

We performed Wolff simulations of the 2D Ising 
model at criticality, using SR with feed-back positions 
(p,g)=(36,ll), (89,38), (127,64) and (250,103) as listed 
in Ref. I51J and references therein. For each pair (p,L) 
we took 100 samples of 10 6 Wolff clusters. Thus we de- 
termined the coefficient in Eq. (ph: p — 1.09(1) Lmax- 
Here, and below, the numbers in parentheses indicate 
the statistical errors. 

The results for the energy deviations 8E = (E/E cx — 1) 
are plotted in Fig. 1. The exact results are taken from 
Ref. The maximum deviations occur at L = 7, 12, 
15 and 22 respectively, in agreement with Eq. (||). The 
inset in Fig. 1 displays the maximum deviations of the 
energy SE max as a function of the shift-register length. 
A fit yields SE max cx p -°-^). 

The resulting data collapse for the scaled deviations 
SE = p Q 88 SE is shown in Fig. 2 versus the scaled system 
size L = p~ 0A3 ( 5 >L. The linear decay on the right obeys 
^ocZ- a84 ( 4 ). 

If the data for L > p 4 / 7 keep following the linear trend 
in Fig. 2, the maximum possible deviations can be de- 
scribed by relation 



5E<0.3L-° M p-° 52 . 



(3) 
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The results for (127,64) do not fit the curve well. This 
is no surprise because shift registers with (p,q) close to 
powers of 2 are known Jl?]] to produce relatively poor 
random numbers. 

Similarly, we sampled the deviation of the specific heat 
C. Fig. 3 shows scaled deviations SC = p 051 ^SC ver- 
sus the scaled system size which is the same as for 8E. 
For large L this curve behaves as SC oc L-°- 21< - 2 \ The 
deviations satisfy 

-^<0.85L- a21 p- a42 (4) 

but they can also be decribed in terms of a logarithm of 
L plus a constant. 

Fig. 4 shows analogous results for the dimensionlcss 
ratio Q = (m 2 ) 2 /(m 4 ), which is related to the Binder 
cumulant |Q, using SQ = SQ p a60 W along the vertical 
scale. On the right hand side the data behave as SQ oc 
£-0.45(5)^ Extrapolation leads to 

SQ< 0.244 L-° A5 p-° A1 . (5) 

In order to explain the origin of the observed devia- 
tions, we present a simple model that captures the essen- 
tials of the Wolff cluster formation process. This model 
simulates a directed random walk in one dimension [ p"9| . 
At discrete times, the walker makes a step to the right 
with probability fi; otherwise the walk ends. The proba- 
bility to visit precisely n consecutive nodes is 

-Pcx(n) = M™ 1 (1 — M)- (6) 

Now, we simulate this model using a SR random-number 
generator. Each walk starts directly after completion of 
the preceding one, without skipping any random num- 
bers. First, we use the 'positive' condition x n > /i for 
stopping. Thus, the random number at start always ful- 
fills the condition xq > ^i, which ended the preceding 
walk. 

In the simplest case [i = 1/2, only the leading bit af- 
fects this condition. As long as the walk proceeds, the 
leading bits of the random numbers x n are zero. After 
p — 1 successful moves, the SR algorithm will produce 
a number x p with the leading bit equal to 1. Thus the 
walker cannot visit more than p nodes. 

A probabilistically equivalent condition for stopping is 
the 'negative' condition x n < 1 — /a. Then, the leading 
bit of xq must be 0, and for x n (n > 1) it is 1 until the 
walk ends. The walk cannot stop at the n = p, since 

Xp £B Xp—q = (£) 1 = 1. 

One can calculate the deviation from the exact value of 
P(n) at n = p, n = p + q and at all linear combinations of 
numbers p and q. The detailed analysis will appear else- 
where and here we only mention that the probability 
deviation 8P{n) = {P comp {n) - P cx (n)) / P cx (n) at n = p 
for the posititive condition is equal to (1 — fx)/ fx. It is im- 
portant that a deviation, even at only one point n = p, 



results in a deviation of the probability function for the 
points n > p by SP(n) = (2/i - l) 2 /^ 4 - 1 < 0. The de- 
viations at the 'resonances' n — ip + jq (i — 1, 2, ... and 
—i<j< i) are positive and lead to negative deviations 
of the next points. 

Thus, in the case of the positive condition, most of the 
8P(n) are negative. In the case of the negative condition, 
5P(n) is negative for n = p 2 k (k = 0, 1, 2, ...); this results 
in positive deviations for the following points. 

In effect, this replaces the probability p by a new 'ef- 
fective' probability fj,*, with fx* > p for the positive con- 
dition and p* < p for the negative condition for most 
n > p. This provides a qualitative explanation of the de- 
viations in Wolff simulations. The completion of a Wolff 
cluster is strongly correlated with the value of the random 
numbers used at that time. Thus, the three-bit correla- 
tions generated by the production rule lead to two-bit 
correlations in the following p random numbers. In par- 
ticular when the mean Wolff cluster size is about p, one 
may expect serious deviations in the calculated quanti- 
ties. 

When one replaces the positive by the negative condi- 
tion, in effect the three-bit correlation is inverted. Thus, 
one expects a change of sign of the systematic errors. We 
confirmed this for the 2D Ising model. 

A simple modification of the SR (1) is to use only 
one out of every m random numbers generated by the 
production rule If m = 2 k , k = (1,2,...) this 

will lead to the same production rule (1). For m = 3 
and, as an example, for SR (36,11) the resulting produc- 
tion rule is (36,24,12,11): a 5-point production rule, i.e. 
x n = x n -3e®x n -24@x n -i2®x n -xx- However, the lowest- 
order correlations of the resulting random numbers do 
not occur at n = p = 36, but at n — 48 because the 
production rule is equivalent with a 4-point one, namely 
(48,23,11) (^o|. The effect due to 4-point correlations ap- 
pears to dominate over the 5-bit effects for p > 1/2. The 
deviations 5P(n) of Eq. (6) resemble those for a 3-point 
production rule. But for /i close to 1 they stand out only 
at n = 48 k, k = 1, 2, .., and not at linear combinations 
of other magic numbers. Their sign is the same for the 
positive and negative conditions because the 4-point rule 
correlates an even number of bits). 

Next, we investigated these 4-bit effects in the case 
of Wolff simulations, using every third number produced 
by the rules (36,11) and (89,38), and runs of 10 9 clus- 
ters. The deviations obey the same scaling laws, but the 
amplitudes are about 20 times smaller for each of the 
quantities E, C and Q, in accordance with the behavior 
of the ID model (see the asterisks in Fig. 2). 

For m = 5 - using only every 5-th number |Tl"| ] - the 
effective production rule correlates 5 bits [^0| . It leads to 
deviations in ID model, in particular at n = pk, k = 2\ 
They are less than for the SR of Eq. (1) Q. 

Very long simulations, using 100 samples x 10 7 Wolff 
steps for m = 5, show that the deviations are even smaller 
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than for m = 3. Table Q displays data for SR (36,11) and 
(89,38) at lattice sizes L = 7 and L = 12, respectively, 
Similar data are included for m = 3 and for m = 1. 

So, we propose, in addition, that the systematic de- 
viations of 2D Ising Wolff simulations are described by 
Eqs. (3-5) for all SR-type algorithms, but the coefficients 
should be corrected with a factor of roughly 10~( mc ~ 3 ), 
where m c is the number of bits correlated by the produc- 
tion rule. 

A preliminary analysis pl| confirms relation (2) also 
for the 3D Ising model. The deviations can also be col- 
lapsed on universal curves, but the exponents and ampli- 
tudes differ from the 2D case. 

We conclude that the ID model provides a useful way 
for the analysis of random numbers, in particular for the 
detection of harmful correlations in SR sequences. The 
errors in Wolff simulations induced by these correlations 
satisfy scaling relations which have a considerable sig- 
nificance for large-scale Wolff simulations. For instance, 
they confirm that in recent simulations [fbif of the ran- 
dom bond Ising model with lattice sizes L greater than 
128, the bias due to the (250,103) Kirkpatrick-Stoll rule 
was less than the statistical errors. 

As explained above, 3-bit correlations in a SR produc- 
tion rule lead to 2-bit correlations in the first p random 
numbers used for the construction of a new Wolff cluster. 
If the size of the latter grows large in comparison with p, 
the 2-bit effect will decrease because the amount of cor- 
relation contained in the first p numbers remains finite. 
Indeed, this is in agreement with the power-law decay on 
the right-hand sides of Figs. 2-4. Although 3-bit effects 
seem to be much smaller in the cases (L,p) investigated 
by us, there is no reason to believe that they are absent. 
Thus, eventually they are expected to end the aforemen- 
tioned power-law decay. 
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TABLE I. Deviations of energy SE, specific heat 8C and 
ratio SQ. The statistical error in the last decimal place is 
shown between parentheses. We used a shift-register length 
p = 36 for L — 7 and p = 89 for L = 12. The bias appears to 
depend strongly on the number m c of bits correlated by the 
production rule. 
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m c 
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SC 


SQ 


7 
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-0.094307 ( 52) 


0.014442 ( 10) 
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-0.000356 ( 13) 


0.005894 ( 69) 
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-0.000060 ( 11) 


0.001122 ( 60) 


-0.000133 ( 15) 


12 
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0.003345 ( 9) 


-0.066797 ( 65) 


0.009577 ( 13) 


12 
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-0.000149 ( 15) 


0.003296 ( 79) 


-0.000274 ( 18) 


12 


5 


-0.000003 ( 11) 


0.000136 ( 89) 


-0.000009 ( 15) 
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FIG. 1. Energy deviations SE for several SR, namely FIG. 3. Scaled deviation of specific heat 8C versus the 

(36,11): o; (89,38): +; (127,64): □; and (250,103): A. The scaled system size, for several SR. The symbols are defined 
inset shows the maximum value of 8E as a function of p. in the caption to Fig. 1. 




FIG. 2. Scaled deviation of the energy SE versus the scaled FIG. 4. Scaled deviation of dimcnsionless ratio 8Q versus 

system size, for several SR. The symbols are defined in the the scaled system size, for several SR. The symbols are defined 
caption to Fig. 1. in the caption to Fig. 1. 
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FIG. 1. Energy deviations SE for several SR, namely (36,11): o; (89,38): +; (127,64): □; and (250,103): A. The inset shows 
the maximum value of SE as a function of p. 
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FIG. 2. Scaled deviation of the energy SE versus the scaled system size, for several SR. The symbols are defined in the 
caption to Fig. 1. 
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FIG. 3. Scaled deviation of specific heat SC versus the scaled system size, for several SR. The symbols are defined in the 
caption to Fig. f . 
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FIG. 4. Scaled deviation of dimensionless ratio SQ versus the scaled system size, for several SR. The symbols are defined in 
the caption to Fig. 1. 
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